subroutine div_cleaning
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use corgan_com_M
      use cindex_com_M
      use numpar_com_M
      use cprplt_com_M
      use cophys_com_M
      use ctemp_com_M
      use blcom_com_M
      use objects_com_M
      use controllo_M
      use celeste3d
      implicit none
!-----------------------------------------------
!   G l o b a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------

!     calculate the source for Poisson's equation
!     the call computes a correction to E0
!
               call divc (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x&
                  , c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y&
                  , c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vol, ex0, ey0, ez0&
                  , divpix)
!
	dnorm=0.
               do n = 1, ncells
                  ijk = ijkcell(n)
                  wate(ijk,10) = fourpi*qdnc0(ijk) - divpix(ijk)
		  dnorm=dnorm+(fourpi*qdnc0(ijk) - divpix(ijk))**2
               enddo
	       write(*,*)'prima=',dnorm
!
!     with these parameters, poisson solves poisson's equation
!
               dumdt = 0.0
               dumscal = 1.
               wk = 0.
               phibc = 0.
!
!
               call Poisson(ncells, ijkcell, nvtxkm, nvtx, ijkvtx, itdim, iwid&
                  , jwid, kwid, precon, tiny, relax, ibp1, jbp1, kbp1, cdlt, &
                  sdlt, strait, dx, dy, dz, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y&
                  , c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z&
                  , c6z, c7z, c8z, vol, vvol, wate(1,28), wate(1,1), wate(1,2)&
                  , wate(1,7), wate(1,8), itmax, error, wate(1,10), dumdt, &
                  dumscal, wk, phibc, ex, ey, ez, wate(1,4), wate(1,5), wate(1,&
                  6), phipot, wate(1,9), wkl, wkr, wkf, wke, periodicx)
!
!
               call gradf (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, &
                  c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, &
                  c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, phipot, gradxb, &
                  gradyb, gradzb)
!
!
               call axisgrad (ibp1, jbp1, iwid, jwid, kwid, gradxb, gradyb, &
                  gradzb) !....nothing is done with this call
!
!dir$ ivdep
               do n = 1, nvtx
                  ijk = ijkvtx(n)
                  rvvol = 1./vvol(ijk)
                  ex0(ijk) = ex0(ijk) - gradxb(ijk)*rvvol
                  ey0(ijk) = ey0(ijk) - gradyb(ijk)*rvvol
                  ez0(ijk) = ez0(ijk) - gradzb(ijk)*rvvol
               enddo

         call divc (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x&
                  , c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y&
                  , c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vol, ex0, ey0, ez0&
                  , divpix)
!
			   dnorm=0.
               do n = 1, ncells
                  ijk = ijkcell(n)
                  wate(ijk,10) = fourpi*qdnc0(ijk) - divpix(ijk)
				  dnorm=dnorm+(fourpi*qdnc0(ijk) - divpix(ijk))**2
               enddo
	       write(*,*)'dopo=',dnorm



end subroutine div_cleaning

